Structural transitions and arrest of domain growth in sheared binary immiscible fluids 

and microemulsions 
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We investigate spinodal decomposition and structuring effects in binary immiscible and ternary 
amphiphilic fluid mixtures under shear by means of three dimensional lattice Boltzmann simulations. 
We show that the growth of individual fluid domains can be arrested by adding surfactant to the 
system, thus forming a bicontinuous microemulsion. We demonstrate that the maximum domain 
size and the time of arrest depend linearly on the concentration of amphiphile molecules. In addition, 
we find that for a well defined threshold value of amphiphile concentration, the maximum domain 
size and time of complete arrest do not change. For systems under constant and oscillatory shear 
we analyze domain growth rates in directions parallel and perpendicular to the applied shear. We 
find a structural transition from a sponge to a lamellar phase by applying a constant shear and 
the occurrence of tubular structures under oscillatory shear. The size of the resulting lamellae and 
tubes depends strongly on the amphiphile concentration, shear rate and shear frequency. 

PACS numbers: 64.75.+g,47.11.-j, 82.70.Uv, 83.60.-a 

I. INTRODUCTION 

Complex fluid mixtures consisting of immiscible fluid 
species and surfactants are ubiquitous in many industrial 
applications. These fluid mixtures involve both hydro- 
dynamic flow effects and complex interactions between 
fluid particles. Typical examples can be found in the 
food, cosmetic and chemical industries where surfactants 
are utilized to stabilize otherwise immiscible fluids. A 
good example is a barbecue sauce containing large frac- 
tions of water and oil or fat. Without any additives 
the constituents would phase separate, entering a not 
very appealing de-mixed state in which the fat accumu- 
lates in a thick layer on top of the remainder. Adding 
an emulsifier or surfactant helps to stabilize the sauce. 
These molecules are often called amphiphiles and, in their 
simplest form, are comprised of a hydrophilic (water- 
loving) head group and a hydrophobic (oil-loving) tail. 
The surfactant molecules self-assemble on the surface of 
oil droplets and reduce the surface tension. Thus, the 
droplets stabilize and remain suspended within the bulk 
water. A typical emulsifier used by the food industry 
is egg yolk lecithin. Proteins and emulsifiers with low 
molecular weight are also common emulsifiers. 

Amphiphilic fluids containing at least one species made 
of surfactant molecules are in general complex fluids 
that can self-assemble to form ordered structures such 
as lamellae, micelles, not ordered sponge phases and liq- 
uid crystalline (cubic) phases. In general, adding am- 
phiphiles to a binary mixture of otherwise immiscible flu- 
ids (for example oil and water) which is undergoing phase 
separation can cause the de-mixing process to slow down. 
If the amphiphile concentration is sufficiently high, the 
de-mixing process might eventually arrest completely. It 
has been shown by Langevin, molecular dynamics, lattice 
gas, and lattice Boltzmann simulations that the temporal 



growth law for the size of oil and water domains in a sys- 
tem without amphiphiles follows a power law t a [l|, and 
crosses over to a logarithmic growth law (In t)® .where a, 
6 are fitting parameters and t is the time 0, EL Q • A fur- 
ther increase of the surfactant concentration can lead to 
growth which is well described by a stretched exponen- 
tial form A — B exp(— Ct D ), where capital letters denote 
fitting parameters [J, Q . By adjusting temperature, fluid 
composition or pressure, amphiphiles can self-assemble 
and force the fluid mixture into a number of equilib- 
rium structures or mesophases. These include lamellae 
and hexagonally packed cylinders, micellar, primitive, di- 
amond, or gyroid cubic mesophases as well as sponge 
phases. In this paper we focus on the sponge mesophase, 
which in the context of our simulations is called a bi- 
continuous microemulsion since it is formed by the am- 
phiphilic stabilization of a phase-separating binary mix- 
ture, where the immiscible fluid constituents occur in 
equal proportions. Here, the oil and water phases in- 
terpenetrate and percolate and are separated by a mono- 
layer of surfactant at the interface. 

Such complex fluids are often subject to shear forces 
and show pronounced rheological properties 0, Q • Of- 
ten, shear induced transitions from isotropic to lamel- 
lar phases can be observed. These have been studied 
experimentally in binary Q and ternary amphiphilic flu- 
ids [{J [l(| • If oscillatory shear is applied, a further transi- 
tion to a tubular phase or a transition between differently 
oriented lamellar phases can occur [ll|, [l2|, 0, [H| • 

Computationally, such complex fluids are too large 
and expensive to tackle with atomistic methods such as 
molecular dynamics, yet they require too much molecu- 
lar detail for continuum Navier-Stokes approaches. Al- 
gorithms which work at an intermediate or "mesoscale" 
level of description have been developed during the last 
twenty years, including dissipative particle dynamics [lH 
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lis |17|, lattice gas cellular automata[18||. the stochastic 
rotation dynamics [Hi , [20L [2l| , and the lattice Boltzmann 
equation [23. I23I l24j| . In particular, the lattice Boltzmann 
method has been found highly useful for simulating com- 
plex fluid flows in a wide variety of systems. This al- 
gorithm, described in more detail in the next section, is 
extremely well suited to implementation on parallel com- 
puters, which permits very large systems to be simulated, 
reaching hitherto inaccessible physical regimes. 

In this paper we investigate spinodal decomposition 
and structuring effects in binary immiscible and ternary 
amphiphilic fluid mixtures under shear by means of large 
scale three dimensional lattice Boltzmann (LB) simu- 
lations. The purely kinetic LB method we use is able 
to model complex flows whose rheological properties are 
emergent from the mesoscopic kinetic processes without 
any imposed macroscopic constraints 25:]. 

Varieties of the lattice Boltzmann method have been 
used successfully to study the behaviour of multi-phase 
flows in the past. A number of authors have investigated 
spinodal decomposition @, HE H3, [H, H, [H HH M, HI, 
Hi], [H, H(| and the same phenomenon has also attracted 
some interest in the presence of shear, where structural 
transitions from isotropic to lamellar or tubular phases 
may occur [3 HI lM \M ■ 

There have been only limited investigations of the in- 
fluence of amphiphiles on the domain growth within the 
lattice Boltzmann method, despite the fact that ternary 
amphiphilic fluids have been studied by a number of au- 
thors [24], [41], [42|, H^|. For example, it has been shown 
that the lattice Boltzmann method can be used to de- 
scribe the self-assembly and the rheological properties of 
mesophases including the primitive P-phase [44j and the 
gyroid phase df| . The gyroid mesophase in particular has 
been of major interest during the last few years, where 
the phase formation and structural properties 0, [46| , the 
influence of defects I471. I481 149| . l50j , as well as its proper- 
ties under shear [45l l5l| have been investigated. 

In this paper, our aim is to focus on the effect of sur- 
factant concentration on the length and time scales of 
arrested growth and on the changes in structural prop- 
erty induced by steady or oscillatory shear. 

The remainder of the paper is organized as follows. 
After an introduction to the simulation method and the 
fluid parameters used, we present results from the simula- 
tion of ternary amphiphilic systems with varying surfac- 
tant density. Here, we draw on the results of Gonzalez- 
Segredo et al. f{| for the simulation parameters used 
within the present work. Thereafter, we extend our study 
to systems under constant and oscillatory shear, where 
we report on the influence of the domain growth rates 
and structural transitions induced by the shear. Finally, 
a summary and conclusions are given. 



II. SIMULATION METHOD 

A set of equations can be used to represent a standard 
LB system involving multiple species [52j 



n?(x + c i ,t + l)-nf(x > t)=n? 



(1) 



with i = 0, 1, . . . , b. The single-particle distribution func- 
tion nf(x,t) indicates the density of species a, having 
velocity Cj, at site x on a D-dimensional lattice of coor- 
dination number b, at timestep t. The collision operator 
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(2) 



represents the change in the single-particle distribution 
function due to the collisions. A popular form is the 
single relaxation time r a , linear 'BGK' form [53[ for 
the collision operator. It can be shown for low Mach 
numbers that the LB equations correspond to a solu- 
tion of the Navier-Stokes equation for isothermal, quasi- 
incompressible fluid flow. The lattice Boltzmann method 
is an excellent candidate to exploit the possibilities of 
parallel computers, as the dynamics at a lattice site re- 
quires only information about quantities at nearest neigh- 
bour lattice sites [24], |4i| . The local equilibrium distri- 
bution nf eq plays a fundamental role in the dynamics of 
the system as shown by Eq. JT]). In this study, we use a 
purely kinetic approach, for which the local equilibrium 
distribution ri* eq (x, t) is derived by imposing certain re- 
strictions on the microscopic processes, such as explicit 
mass and total momentum conservation 1541 
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where p a (x, t) = J2i V? ( x ? t) is the fluid density and 
u = u(x, t) is the macroscopic bulk velocity of the fluid, 
given by p Q (x, t)u a = J2i n i( x ^) c i- & are the coeffi- 
cients resulting from the velocity space discretization and 
c s is the speed of sound, both of which are determined 
by the choice of the lattice. We use a D3Q19 imple- 
mentation, i.e., a three dimensional lattice with 19 dis- 
crete velocities. Immiscibility of species a is introduced 
in the model following Shan and Chen [55l |56| |. where 
only nearest neighbour interactions among the species 
are considered. These interactions are described by a 
self-consistently generated mean field body force 

F a (x, t) ee -^ Q (x, t) J2 g a& E ^( x '> f )( x ' - x ) . ( 4 ) 



where ip a (yL,t) is the so-called effective mass, which can 
have a general form for modeling various types of fluids 



(we use ip a = (I - 



)[55j), and g aB , is a force coupling 



constant whose magnitude controls the strength of the 
interaction between components a, a and is set positive 
to mimic repulsion. The dynamical effect of the force is 
realized in the BGK collision operator by adding to the 
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velocity u in the equilibrium distribution of eq. ([3]) the 
increment 
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(5) 



Amphiphilcs are introduced within the model as de- 
scribed in [25| and [43| . An amphiphile usually possesses 
two different fragments, each having an affinity for one 
of the two immiscible components. The orientation of 
any amphiphile present at a lattice site x is represented 
by an average dipole vector d(x, i). Its direction is al- 
lowed to vary continuously and no information is spec- 
ified for each velocity Cj, for reasons of computational 
efficiency and simplicity. The amphiphile density at a 
given site is given by an additional fluid species with 
density p s . The model has been used successfully to 
study spinodal decompositio n 12. [271 . the formation of 
mesophases 0, 0, Ell IH, HI, l5lll57j. and flow in porous 
media |39j |. 

We use LB3D[49|, a highly scalable parallel LB code, to 
implement the model. The very good scaling of our code 
permits us to run all our simulations on multiprocessor 
machines and computational grids in order to reduce the 
length of data collection to a few weeks. Also, we are 
able to use simulation boxes typically eight times bigger 
than previous studies so as to minimize finite size effects. 

In order to study the rheological behaviour of com- 
plex fluid mixtures, we have implemented Lees-Edwards 
boundary conditions, which were originally developed for 
molecular dynamics simulations [58(. They reduce finite 
size effects as compared to moving solid walls [58| and 
have been used in lattice Boltzmann simulations by var- 
ious authors [37], [H, H3, . This computationally con- 
venient method imposes new positions and velocities on 
particles leaving the simulation box in the direction per- 
pendicular to the imposed shear strain while leaving the 
other coordinates unchanged. Choosing z as the direction 
of shear and x as the direction of the velocity gradient, 
we have 



(z + A z ) mod N z , x > N x 
z mod N z , < x < N x 

(z - A z ) mod N z , x < 



+ U ,x> N. 



- U 



0<x<N x 
x<0 



(6) 
(7) 



where A z = U At, U is the shear velocity, u z is the 
z— component of u and N x r z \ is the system length in the 
x{z) direction. We also use an interpolation scheme sug- 
gested by Wagner and Pagonabarraga [38| as A z is not 
generally a multiple of the nearest neighbour lattice dis- 
tance. For oscillatory shear, we set 



U(t) = Ucos(ojt) 



(8) 



where lu/2tt is the frequency of oscillation. 

In non-sheared studies of spinodal decomposition it has 
been shown that large lattices are needed to overcome 



finite size effects. There, I28 3 was the minimum ac- 
ceptable number of lattice sites Q . More quantitatively, 
Kendon et al.[33[ set L/4 as the maximum length scale 
which is not affected by finite size effects in their spinodal 
decomposition simulation, where L is the length of the 
simulation box. We therefore choose 256 3 for all non- 
sheared simulations to limit the influence of finite size 
effects even further. For high shear rates, systems also 
have to be very extended in the direction of the applied 
shear because, if the system is too small, the domains 
interconnect across the z = and z = N z boundaries to 
form continuous lamellae in the direction of shear [39l.l49l]. 
Such artefacts need to be eliminated from our simula- 
tions. In this case, a good compromise to limit finite 
size effects and to keep the computational expense as 
low as possible is a lattice size of 128x128x512 and this 
is used here. Mass and relaxation times are always set to 
unity, i.e. r"=I.O, to" =1.0. We call the two immiscible 
fluids "red" and "blue" and set their initial densities to 
identical values, p r = p b . The initial average surfactant 
density p s is varied between 0.0 and 0.7. The lattice is 
than randomly populated with constant initial total fluid 
densities p tot = p r + p b + p s = 1.6. This is in contrast to 
previous studies where only p r + p b was kept constant 
The coupling constant in Eq.H] between "red" and "blue" 
species is set to gbr=0.08, the coupling between an immis- 
cible fluid and surfactant to gb s =-0.006 and the constant 
describing the strength of the surfactant-surfactant in- 
teraction is kept at g ss — —0.003. All units in this paper 
are given in lattice units if not stated otherwise. While 
our method has been used to simulate other mesophases 
like lamellar phases, the primitive P-phase (44j, and the 
gyroid phase [f| |4g|, the parameters used in all simu- 
lations presented here are known to produce a sponge 
phase in the absence of bulk flow. More detailed inves- 
tigations of the particular choice of coupling constants 
and how they modify the system's properties are given 

in @,i \M,m mi. 

To analyze the behaviour of the various simulations, 
we define the time dependent lateral domain size L(t) 
along direction i = x, y, z as 



where 



2tt 



£ k S(k,t) 



(9) 



(10) 



is the second order moment of the three-dimensional 
structure function 



5(k,t)^i|^(t)| 2 



(11) 



with respect to the Cartesian component i, () denotes 
the average in Fourier space, weighted by S(k, t) and V 
is the number of nodes of the lattice, (j>' k (t) the Fourier 
transform of the fluctuations of the order parameter 



4 



cj>' = <j>— {(f>), and ki is the ith component of the wave vec- 
tor. A projection of the structure function allows us to 
compare simulation data to scattering patterns obtained 
in experiments. We obtain those projections by sum- 
ming up <S(k, t) in one of the Cartesian directions. For 
example, for the projection in the ^-direction this leads 
to S z (k x ,k y ,t) = J2k 2 S(k,t). 

III. RESULTS 
A. Ternary amphiphilic systems without shear 

Spinodal decomposition of a binary immiscible fluid 
mixture has been studied extensively within our model by 
Gonzalez-Segredo et al. Q. The authors report domain 
size scaling as expected for a crossover from diffusive be- 
haviour to hydrodynamic viscous growth, i.e. the domain 
size grows as I a i 7 , with 7 being between 1/3 and 1. 
Moreover, they find very good agreement with the dy- 
namical scaling hypothesis, recovering the expected uni- 
versal behaviour of the structure function. 

If one adds surfactant to a binary immiscible fluid mix- 
ture, the surfactant molecules self-assemble at the inter- 
face between the two species and slow down the phase 
separation process. For sufficiently high surfactant con- 
centrations, domain growth is arrested completely lead- 
ing to a stable microemulsion. In [5], Gonzalez-Segredo 
et al. extend their work to ternary amphiphilic fluids 
and study how the phase separation of a binary immis- 
cible fluid mixture is altered by the presence of surfac- 
tant. As already described in the introduction, by grad- 
ually increasing the surfactant concentration they find 
a slowing down of the domain growth, initially from 
algebraic to logarithmic temporal dependence, and, at 
higher surfactant concentrations, from logarithmic to 
stretched-exponential behaviour. They also observe a 
structural transition from sponge to gyroid phases by 
increasing the amphiphile concentration or varying the 
amphiphilc-amphiphile or amphiphile-binary fluid cou- 
pling constants. For growth-arrested mesophases, they 
observe temporal oscillations of the domain size due to 
Marangoni flows. 

In the present work we use simulation parameters 
which differ from previous studies and which are known 
to produce a sponge phase. In order to avoid effects due 
to variations of the fluid densities, we also keep the total 
density in the system constant at 1.6 (in lattice units), 
while varying the surfactant densities p s between 0.00 
and 0.70. Furthermore, our simulation lattices are up 
to eight times larger in order to keep the influence of fi- 
nite size effects to a minimum and we simulate for up to 
30000 timesteps in order to gain a better understanding 
of the long time behaviour of the system. We study the 
influence of the amphiphile concentration on the phase 
separation process in this section and reproduce previ- 
ous results with the present parameters. In addition, we 
study the dependence of the maximum achievable domain 



size in a stable microemulsion on surfactant concentra- 
tion as well as the time needed to achieve this state. 

To depict the influence of the surfactant density on the 
phase separation process, Fig. [T] shows three volume ren- 
dered 256 3 systems at surfactant densities 0.0 (left), 0.15 
(center), and 0.3 (right). As in all figures throughout 
the paper, for better visibility only one of the immiscible 
fluid species is shown. Different colours denote the inter- 
face and areas of high density of the rendered fluid. The 
surfactant particles are aligned at the interfaces and the 
second immiscible constituent fills the void space. After 
30000 timesteps the phases have separated to a large ex- 
tent when no surfactant is present (left). Running the 
simulation for even longer would result in two perfectly 
separated phases, each of them contained in a single do- 
main only. If one adds some surfactant (p s — 0.15, cen- 
ter), the domains grow more slowly, visualized by the 
smaller structures in the volume rendered image. For 
sufficiently high amphiphile concentrations (p s = 0.30, 
right) the growth process arrests leading to a stable bi- 
continuous microemulsion with small individual domains 
formed by the two immiscible fluids. 

The projected structure function S z {k x ,k y ,t) ("scat- 
tering pattern") is given in Fig. [2] for two surfactant den- 
sities p s =0.00 (a) and 0.30 (b) at timestep t = 10000. 
As can be clearly seen in Fig. [2^), a strong peak oc- 
curs for small values of k x , k y depicting the occurrence of 
length scales which are of the order of the system size. 
For p s =0.3Q, however, the peaks are by a factor of 100 
smaller and shifted to larger values of k x ,k y . We find a 
volcano-like scattering pattern indicating the dominance 
of small length scales. Since our system is cubic and no 
shear is applied, the projections of the structure function 
in x and y direction are equivalent. 

To investigate the influence of surfactant more quan- 
titatively, in Fig. [3] the time dependent lateral domain 
size L(t) as given in Eq. [5] is shown for a number of sur- 
factant densities p s between 0.00 and 0.50. Since the 
lattice is cubic here, L(t) behaves identically in all three 
directions. Fig. and (b) show identical data, but 
different scalings of the axes. In (a), we plot the data 
linearly in order to give a better impression of the time 
dependence of the growth dynamics. However, in order 
to check which data is best fitted by the various growth 
laws, we provide a log-log scale plot of the same data 
in Fig. E^b). For the first few hundred timesteps, the 
randomly distributed fluid densities of the initial system 
configuration cause a spontaneous formation of small do- 
mains resulting in a steep increase of L(t). For p s =0.00 
domain growth does not come to an end until the do- 
mains span the full system. By adding surfactant we can 
slow down the growth process and for high surfactant 
densities p s > 0.25, the domain growth stops after a few 
thousand simulation timesteps. By adding even more 
surfactant, the final average domain size becomes very 
small and does not grow beyond 7.7 lattice sites. We fit 
our numerical data with the corresponding growth laws 
and find that for p s smaller than 0.15 L(t) is best fit by a 
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FIG. 1: (Color online) Volume rendered fluid densities of 256 3 systems at t=30000 for surfactant densities p s =0.00, 0.15, 0.30 
(from left to right) . For better visibility only one of the immiscible fluid species is shown. Different colours denote the interface 
and areas of high density of the visualized fluid. The surfactant particles (not shown) are aligned at the interfaces and the 
second immiscible fluid component fills the void space. After 30000 timesteps the phases have separated to a large extent if no 
surfactant is present (left). Adding a small amount of surfactant (p s = 0.15, center) causes the domains to grow more slowly, 
as depicted by the smaller structures in the volume rendered image. For sufficiently high amphiphile concentrations (p s = 0.30, 
right) the growth process arrests with the formation of a stable bicontinuous microemulsion. 
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FIG. 2: (Color online) Projected structure function 
S z (k x ,k y ,t) for (a) p s =0.00 and (b) 0.30 at timestep t = 
10000. For the case without surfactant, a strong peak occurs 
for small values of k x ,k y reflecting the dominance of length 
scales which are of the order of the system size. For p s =0.30, 
only much smaller peaks occur for larger values of k x , k y indi- 
cating that only small length scales are present. All quantities 
are expressed in lattice units. 



function proportional to t a . For p s being 0.15 or 0.20, a 
logarithmic behaviour proportional to (\nt) e is observed. 
Increasing p s further results in L(t) being best described 
by a stretched exponential. These results correspond well 
with the findings in [f|. 

We study the dependence of the final domain size 
Lmax(p s ) on the amount of surfactant as depicted in 
Fig. IDJa). It can be observed that the maximum domain 
size decreases linearly from 20.9 for p s = 0.25 with in- 
creasing p s until a threshold value is reached at p s = 0.5, 
where L max (t) = 7.7. Then, L max (p s ) decreases much 
more slowly and stays almost constant. The slope of 
the linear regime corresponds to -52.8. The behaviour of 
imax(pf) and t a rrest(yO s ) is consistent with previous lattice 
gas [J, [60| and lattice Boltzmann studies [j| , where the 
authors determine the dependence of the surface tension 
at a planar interface between two immiscible fluid species 
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FIG. 3: (Color online) Average domain size L(t) for various 
surfactant densities p 3 =0.00, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 
0.40, 0.50. Axes are linear in (a) and logarithmic in (b). 
After the initial spontaneous formation of domains, domain 
growth can be described by a power law. With increasing 
p s domain growth slows down and eventually comes to an 
end at a maximum domain size L m ax(p s ). All quantities are 
expressed in lattice units. 



on the surfactant concentration. In 5], the authors find a 
linear dependence between surface tension and surfactant 
density, but they did not study such high concentrations 
as in the current study. In [601 ] . the surface tension ap- 
proaches zero for high concentrations, i.e. a saturation 
occurs causing the size of the individual fluid domains to 
saturate as well. These effects can be explained as fol- 
lows: adding surfactant to a binary fluid mixture causes 
the amphiphilcs to minimize the free energy in the system 
by assembling at the interface between the two immisci- 
ble fluid species. An increase of surfactant concentration 
causes the interfacial area to be maximised in order to 
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FIG. 4: (a) Maximum domain size L lnax (p s ) and (b) time of 
arrest t a rrest(p s ) for various surfactant densities p s =0.25, 0.30, 
0.35, 0.40, 0.50, 0.7. I/ m ax(p") as well as t^ Trest (p s ) decrease 
linearly with the surfactant density p s and saturate for p s > 
0.5. All quantities are expressed in lattice units. 



accommodate as much surfactant as possible. The in- 
creasing interfacial area causes the individual domains 
to become smaller and i max (p s ) decreases. If the surfac- 
tant concentration becomes very high (p s > 0.5 in our 
case), L max (p s ) saturates due to the maximum possible 
interfacial area being reached and all available area being 
covered with surfactant molecules. More amphiphilcs ac- 
cumulating at the interface would lead to very steep and 
energetically unfavourable gradients of surfactant den- 
sity in the system. Therefore, further amphiphilcs have 
to reside within the bulk fluid phases forming micellar 
structures. Within our model the minimum final domain 
size corresponds to 7.7 lattice sites. However, this value 
can be varied by tuning the coupling constants for the 
amphiphilc-amphiphilc or amphiphilc-fluid interactions. 
A more thorough investigation of the influence of the 
particular choice of the coupling constants on the final 
domain size is of particular interest and shall be investi- 
gated within a future project. 

In Fig. [Hb), the number of simulation timesteps 
£arrcst(/0 s ) needed to reach the final domain size is plot- 
ted. Since the time it takes for the system to relax to its 
equilibrium state directly depends on the final domain 
size, it is consistent with the data presented in Fig. [Ha) 
that a linear dependence of t alTes t(p s ) on the surfactant 
concentration can be observed. While for p s = 0.25 7000 
timesteps are needed to reach the maximum possible do- 
main size, for p s = 0.5 500 timesteps are sufficient. For 
p s > 0.5, t arrcst (p s ) decreases much more slowly than for 
p s < 0.5. The slope of t arrcst (p s ) in the linear regime is 
given by -26000. 
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FIG. 5: (Color online) Domain size L(p s ) in x (a), y (b), and 
z direction (c) for surfactant densities p s =0.0, 0.1, 0.2, 0.3, 
0.4 and a constant shear rate 7 = 1.56 x 10~ 3 . All quantities 
are reported in lattice units. 



B. Steadily sheared systems 

If a binary immiscible fluid mixture is driven mechan- 
ically by external shear forces, it is known that the evo- 
lution of domains and phase separation processes are 
changed profoundly [j| H, [l(| • The most noticeable effect 
is the formation of a lamellar phase, i.e. elongated do- 
mains or lamellae form and align along the flow direction. 
Due to the anisotropy of the system, the time dependent 
domain size L{t) behaves differently for discerned coor- 
dinate axes in this case. Furthermore, modified growth 
exponents are expected due to the anisotropic effects. 

As already seen in the previous section, adding am- 
phiphilcs to a binary immiscible fluid under shear can 
change its properties dramatically. The amphiphiles sta- 
bilize the interface between the immiscible fluid species 
and the domain growth is hindered as described in the 
previous section. Moreover, the amphiphiles might form 
complex structures which can have an impact on the 
properties of the sheared fluid leading to non-Newtonian 

flow HEIIH3- 

We study ternary 128x128x512 sized systems under 
constant shear. The shear rate is set to 7 = 1.56 x 10~ 3 
and 7 = 3.12x 10~ 3 , while the surfactant density is varied 
between p s =0.0 and 0.4. Fig. [5] shows the time depen- 
dent lateral domain size for all three coordinate axes at 
7 = 1.56 x 10~ 3 . In the x direction, which is the axis 
between the shear planes, the power law regime of L x (t) 
starts at £=500 for the p s =0.0 case, while for higher p s 
the initial growth regime is overcome by the power law 
regime before the first measurement at £=100. As long 
as p s < 0.3, the growth rate is not hindered by the am- 
phiphiles and domains grow until the end of the simula- 
tion. For p s =0.3 the power law regime starts to break 
down at £=900 and L x (t) saturates at £=5000. Adding 
even more surfactant results in an even earlier saturation 
at £=1500. The y direction is the direction parallel to the 
shear planes and perpendicular to the direction of shear. 
Since this direction is less affected by the shear forces, 
L y (t) grows faster than L x (t) for low surfactant concen- 
trations p s <0.2 causing the domains forming being ex- 
tended in the y direction. For p s =0.2 L x (t) and L y (t) 
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FIG. 6: (Color online) Volume rendered (128x128x512 simulation boxes) fluid densities for surfactant density p s =0.2, a constant 
shear rate j = 1.56 x 10 -3 and variable number of time steps £=1000 (upper left), £=4000 (upper right), £=6000 (lower left), and 
£=10000 (lower right). While the shape of individual domains does not show distinct features at early times of the simulation, 
elongated structures appear at £=4000 and start to become aligned with the shear at £=6000. At late stages of the simulation 
run (£=10000), the lamellar phase characterised by thin and long lamellae filling the whole system can be observed. 



behave almost identical, while for p s >0.2 a crossover 
occurs and the maximum attainable value for L y (t) is 
below the result for L x (t). In the direction of shear (z 
direction), L z (t) saturates even for the no surfactant case 
at L z (t)=25 and comes to arrest at even smaller values 
with increasing p s . The complex behaviour of Li(t) can 
be better understood by reminding ourselves that the 
domain size is measured in the direction of the Carte- 
sian coordinate axes. However, individual fluid domains 
occurring in the system are being elongated due to the 
shear and try to align with the shear gradient. Thus, 
they are not parallel to any coordinate axis. Therefore, 
with a measurement of L z (t) one is not able to detect 
the actual length of individual lamellae, but only their 
thickness in the z direction. Similar arguments are valid 
for L x (t), shear causing the measured domain size in the 
x direction to be larger than the lamellae's thickness. For 
increasing p s , the average domain size reduces due to the 
influence of the amphiphiles, thus causing the individual 
domains to become smaller. If p s >0.2, the alignment of 
the domains with the shear causes L x (£) to appear larger 
than L z (t). For high surfactant concentrations (p s =0.4) 
all three directions behave very similarly: domain growth 
comes to an end after less than 2000 timesteps and the 
final domain size is between 10 and 15 lattice units in 
all directions , signaling the appearance of a stable mi- 
croemulsion where the shape of the domains is almost 
unaffected by the shear. 

Regular peaks occur in L z (t) at every 2500 timesteps 
with less pronounced peaks in between them. These 
peaks can be explained as follows: For the stretching 
of domains, a certain amount of work against surface 
tension is needed. On macroscopic scales, the stress 
tensor does not vanish due to the viscoelastic response 



of the system [gj HH. On the microscale, however, 
a breakup and recombination of domains can be ob- 
served [63]. These domains grow by diffusion and even- 
tually join each other to form larger structures. If the in- 
ternal stress becomes too large due to the shear induced 
deformation, they break up and start to form again. As- 
suming a large system with many independent domains 
growing and breaking incoherently, the only observable 
effect might be a slowing down of the domain growth. In 
contrast, if the growth and breakup occur coherently as 
they do in our simulations, a periodicity in the measured 
time dependent domain size can be observed j&ij . As can 
be observed in Fig. [5J:) , the frequency of domain breakup 
is independent of the surfactant concentration, while the 
height of the peaks decreases with increasing p s . 

Figure [6] shows volume rendered examples of a simu- 
lated system with surfactant density p s =0.2 and a con- 
stant shear rate of 7 = 1.56 x 10~ 3 . The four snapshots 
are taken a different times £=1000 (upper left), £=4000 
(upper right), £=6000 (lower left), and £=10000 (lower 
right). It can be observed that at early stages of the sim- 
ulation, the shape of individual domains does not show 
distinct features, while at £=4000, slightly elongated do- 
mains start to occur which begin to align with the shear 
gradient. At £=6000, these features are substantially 
more dominant and at late simulation times (£=10000) 
the system is filled with elongated and thin lamellae con- 
sisting of one of the immiscible fluid species and which 
are almost parallel to the shear plane. 

In order to permit a comparison with experimentally 
available scattering data, in Fig. [7] we present projected 
structure functions for the surfactant-less case (a-c) and 
a surfactant density of p s = 0.30 (d-f) at £ = 10000. 
The x, y and z directions are shown from top to bottom. 
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FIG. 7: (Color online) Projected structure functions ("scat- 
tering pattern") Si for surfactant concentrations p s — 0.0 
(a-c) and p s — 0.3 (d-f) at t — 10000. Projections are in 
x-, y- and z directions (from top to bottom) and for a cubic 
cutout of the full system with side length iV=128. 



The shown projections are for a cubic 128 3 cutout of the 
elongated systems. In contrast to the non-sheared case, 
all three directions show distinguished properties: For 
p s = 0.00, a high peak of S z (k x ,k y ,t) can be observed 
around k x = k y = 0, while in x- and y direction two 
less high peaks at positions > show up. These data 
can be interpreted as follows: At t — 10000, the domain 
size in the direction of the flow corresponds to the size 
of the cubic cutout, i.e. 128 lattice sites. In the x and 
y directions, however, the size of the occurring struc- 
tures is smaller, indicating the occurrence of very long 
lamellar structures in the system. By adding surfactant 
to the system (Fig. EI d-f)), the occurring length scales 
depicted decrease by the splitting of the single peak in 
the z-projection and the occurrence of two small peaks 
at k x = and k y = ±10. While the peaks in the y di- 
rection denote similar length scales as in the z direction, 
the projections in the direction between the shear planes 
(x) show a different behaviour. Here, two parallel struc- 
tures at k y = ±10 and k z between -20 and 20 indicate 
a much wider variation of the thickness of the individual 
domains. This is in contrast to the non-sheared case in 
Fig. [2J where a volcano-like shape of the structure factor 
was observed. 

Doubling the shear rate to 7 = 3.12 x 10 -3 results in 



FIG. 8: (Color online) Domain size L(p s ) in (a) x-, (b) y-, 
and (c) z direction for surfactant densities p 3 — , 0.1, 0.2, 0.3, 
0.4 and a constant shear rate 7 = 3.12 x 10 -3 . 



very similar behaviour, as shown in Fig. [H] In the z di- 
rection, peaks can now be observed even for p s = 0.4, 
but L z (t) is much more noisy for lower surfactant con- 
centrations. However, it can still be seen that there is 
a number of equidistant peaks for p s < 0.4 which oc- 
cur every 2500 timesteps with some additional peaks in 
between in the case of p s = 0.0 and 0.1. The equidis- 
tant peaks occurring with the same frequency as in the 
7 = 1.56 x 10~ 3 case shows that the breakup phenomena 
observed are independent of the shear rate. 

A number of experiments have reported that the shear 
stabilizes the system and causes the phase separation to 
come to an end with LJt) being very large and L x (t) 
being much smaller [3(| Eil [H, [fill ■ We are not able to 
reproduce these results due to the limited size of our sim- 
ulations: substantially larger systems and higher shear 
rates need to be studied in order to quantify the arrest 
of domain growth due to shear. However, the computa- 
tional resources needed would be at the limit of what can 
be done on current supercomputers. 

We have shown in this section that the dynamical scal- 
ing hypothesis does not hold for sheared ternary systems 
in three dimensions since we indeed find three individual 
length scales pointing out the transition from the sponge 
to the lamellar phase: while in the flow direction (2), 
L(t) is determined by the resultant length of the occur- 
ring lamellae, in the direction between the shear planes 
(x) , the domains grow steadily and exhibit power law be- 
haviour up to a maximum that depends on the surfactant 
concentration. In the y direction, domain growth is not 
hindered by shear. In fact, L y (t) grows slightly faster 
than in the non-sheared case. Increasing the surfactant 
concentration has a strong impact on domain growth: 
starting at p s =0.3, L y (t) and L z (t) recover the behaviour 
of the case without shear, i.e. the length scales saturate 
around 15. In the x direction, however, growth contin- 
ues to up to L x (t)=26. This can be explained as fol- 
lows: with increasing surfactant concentration, the final 
domain sizes become smaller, reducing the influence of 
the shear forces in the y and z directions. In the direc- 
tion between the shear planes, however, an increase of 
L x (t) can still be observed because the domains are still 
being elongated due to shear and try to align with the 
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velocity profile. Thus, they are tilted and their size ap- 
pears to be smaller than it actually is in z direction and 
larger in the x direction. 

Our findings are in agreement with Ginzburg-Landau 
and Langevin calculations [H, [6?], [H, H^] as well as 
two-dimensional lattice-Boltzmann simulations of bi nary 
immiscible fluid mixtures as presented in [3l], [H, l40j |. 
However, to the best of our knowledge, there are no de- 
tailed theoretical studies of the dependence of domain 
growth properties on the surfactant concentration. The 
only known work utilizes a Ginzburg-Landau free-energy 
approach to study sheared microemulsions, but does not 
vary the amount of surfactant. In addition, the authors 
only cover two-dimensional systems and are thus unable 
to describe the behaviour of L y (t) [zjj. 



C. Complex fluids under oscillatory shear 

In the case of oscillatory shear, the morphology and 
the domain growth are altered significantly, although the 
average deformation is zero after each period of shear. 
For example, it has been found experimentally for binary 
fluid mixtures that for very low oscillation frequencies 
domain growth can be interrupted [Til ] . or domains can 
grow on much longer time scales than given by the oscil- 
lation frequency [721 ] . Simulations so far either do not in- 
clude hydrodynamic effects, or are limited to two dimen- 
sions [14]. It has been observed in the two-dimensional 
lattice Boltzmann studies by Xu et al. that hydrody- 
namic effects must not be neglected in the case of os- 
cillatory shear since there exists a finite time inversely 
proportional to the viscosity which is required to set a 
linear velocity profile in the system. For high oscillation 
frequencies, this time is longer than the oscillation pe- 
riod and there will never be a linear velocity profile in 
the system, thus influencing the domain growth substan- 
tially [3. 

Due to its higher complexity, oscillatory shear has 
much less been a subject of research projects than sys- 
tems under steady shear. Therefore, the number of pub- 
lications that can be found in the literature is substan- 
tially smaller. To the best of our knowledge, no detailed 
three-dimensional simulation studies of phase separation 
in ternary amphiphilic fluid mixtures under oscillatory 
shear have been reported on so far. The most detailed 
three-dimensional simulations we are aware of are our 
own studies of the gyroid mesophase under oscillatory 
shear as presented in [45| . 

In this section we present our results of simulations 
of systems equivalent to the ones considered in the pre- 
vious section, but with the shear plane moved as given 
by Eq. [51 We apply two different oscillation frequencies 
uj = 0.001 and uj = 0.01, where a single oscillation takes 
6283 timestcps in the slow case and 628.3 timesteps in 
the fast case. 

Let us first consider the case with a lower oscillation 
frequency and lower shear rate, i.e. lu — 0.001 and 



7 = 1.56 x 10 -3 . In the case of oscillatory sheared sys- 
tems, the individual fluid domains try to align with the 
velocity gradient as in the previous section. However, 
since we do not consider steadily moving shear planes 
here, domains are never able to reach a steady state and 
instead have to follow the oscillation of the planes. The 
frequencies considered in our simulations are compara- 
bly high since no linear velocity gradient sustains long 
enough during a single oscillation for the domains to fully 
align with it. This is depicted in Fig. [5] which shows 
two typical examples from a simulation with p"=0.2. On 
the left hand side, a volume rendered snapshot is given 
at £=2500. Here, the oscillating shear planes have just 
passed their reversal point. Close to the shear planes, 
the domains are aligned vertically because they have to 
be turned around in order to follow the changing direc- 
tion of movement of the shear planes. In the bulk of the 
system, however, no preferred direction can be observed 
since the velocity gradient does not interpenetrate the 
whole system. At t= 10000, the shear planes arc in a 
position just before their reversal point. Thus, the fluid 
mixture was accelerated for more than 2000 timesteps 
and the domains close to the shear boundary are well 
aligned in the direction of the flow. In the bulk, again no 
preferred direction can be observed. 

The time dependent lateral domain size of this simula- 
tion and for varied surfactant concentration is presented 
in Fig. |T0j As in the case of continuous shear, the do- 
main growth in y direction is almost uninfluenced by the 
applied shear forces. In fact, for low surfactant concen- 
trations p s < 0.2, L y (t) even grows slightly faster than in 
the case without oscillatory movement. For p s > 0.2 the 
maximum domain size obtained is similar to the steady 
shear case. Due to the non-steady movement of the shear 
planes, L x (t) and L z (t) show a more rich behaviour: both 
functions show distinct kinks around the reversal points 
of the shear and for p s > 0.2 it is found that the growth 
rates are smaller than in the case of steady shear. Thus, 
we can observe the formation of tubular structures which 
are elongated in the y direction and show similar length 
scales in x and z direction. For very high surfactant 
concentrations (p s =0.4) it is not possible to distinguish 
between tubular and spherical structures due to the small 
size of the individual domains. 

In the z direction we can still observe peaks related 
to the formation and breakup, as well as the rotational 
movement, of domains. However, due to the overlaid 
effect of the oscillation, these peaks are not equidistant 
as in the steady shear case anymore. 

In Fig. [TT] we present a system with a ten times larger 
oscillation frequency, i.e. uj — 0.01. Here, the frequency 
of the oscillations is so high that the fluid is not able to 
follow the movement of the walls anymore. Thus, the 
influence of the shear on the growth behaviour becomes 
less pronounced, with the domains constantly growing as 
long as the amount of surfactant present allows it. The 
growth rates are comparable to the non-sheared case here 
and show identical growth laws as in the non-sheared 
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FIG. 9: (Color online) Volume rendered fluid densities for surfactant density p s —0.2 at £=2500 (left) and t=10000 (right). The 
shear rate is 7 = 1.56 x 10 -3 and u> = 0.001. At £=2500, the shear velocity is close to its reversal point and the domains are 
aligned vertically close to the shear plane, while in the bulk of the system no preferred orientation can be observed. After 10000 
timesteps, however, the domains close to the shear plane are aligned with the flow direction because the shear planes are in a 
position just before their reversal point. In the bulk, still no preferred orientation can be found since it takes longer for the 
velocity gradient to penetrate the whole system than the duration of a single period of shear. 
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FIG. 10: (Color online) Domain size L(p s ) in x- (a), y- (b), 
and z direction (c) for surfactant densities p s =, 0.1, 0.2, 0.3, 
0.4 and oscillatory shear with 7 = 1.56 x 10~ 3 , uj = 0.001. 
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FIG. 11: (Color online) Domain size L(p a ) in x (a), y (b), 
and z direction (c) for surfactant densities p 3 =, 0.1, 0.2, 0.3, 
0.4 and oscillatory shear with 7 = 1.56 x 10~ 3 , uj = 0.01. 



case. The only difference is that the exponents are found 
to be smaller while L y (t) grows slightly faster than L x (t) 
and L z (t) depicting the occurrence of tubular structures 
in the system. The z direction is the only component 
of the time dependent lateral domain size that differs 
from the unsheared case, because strong oscillations start 
to appear due to the distortions caused by the moving 
boundaries. 

We have increased the shear rate to 7 = 3.12 x 10~ 3 , 
but the lateral domain sizes obtained are almost identical 
to the 7 = 1.56 x 10~ 3 case. Therefore, we do not present 
an additional figure, but our findings can be used to argue 



that for high oscillation frequencies, the system behaves 
in a manner equivalent to the non-sheared case. 

In the case of oscillatory shear we have shown the 
occurrence of tubular structures due to shear imposed 
anisotropic domain growth, the slowing down of the do- 
main growth rate depending on the oscillation frequency, 
as well as that a microemulsion with high surfactant con- 
centration stays unaffected by external shear forces. Our 
results are in agreement with the simulations of Qiu et 
al. [ll| and the two dimensional lattice Boltzmann sim- 
ulations of Xu et al. 1141. 



IV. CONCLUSIONS 

In this paper we have presented for the first time de- 
tailed three-dimensional lattice Boltzmann studies of bi- 
nary immiscible and ternary amphiphilic fluid mixtures 
under constant and oscillatory shear. 

We have reproduced the well-known power law growth 
of domains in the case of binary immiscible fluids (spin- 
odal decomposition) which crosses over to a logarithmic 
law and to a stretched exponential if one increases the 
surfactant concentration even further. For sufficiently 
high surfactant concentrations, domain growth can come 
to an end and the system corresponds to a stable bi- 
continuous microemulsion. For amphiphile concentra- 
tions of up to 30% we find linear dependencies of the 
time of arrest as well as the maximum domain size on 
the amphiphile concentration. For concentrations above 
30%, both arrest length and time of arrest do not change 
any more since the surface tension at the interfaces be- 
tween the two immiscible fluids is at its minimum. The 
whole interface is filled with amphiphiles and further am- 
phiphile molecules have to reside within the bulk fluid. 

In sheared systems, we have studied the influence of 
moving boundaries on the effect of domain growth and re- 
port domain breakup phenomena depending on the shear 
rate as well as the amphiphile concentration. Depending 
on the surfactant concentration and the shear rate, we 
find a transition from a sponge phase to a lamellar phase. 
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Under oscillatory shear and with the oscillations fre- 
quencies chosen, no linear velocity gradient can build 
up within a single period of shear. Thus, the domains 
are constantly rearranging and align with the flow in the 
vicinity of the shear planes. In the bulk, however, no pre- 
ferred alignment can be observed. But since the growth 
in the y direction is not hindered by the shear, tubular 
structures occur and are best observable for low surfac- 
tant concentrations. For very fast oscillations (u> = 0.01), 
the system is not able to follow the external shear at all. 
Thus, it behaves similar to a non-sheared one. The only 
differences are that fluctuations in L z (t) can be observed 
due to the oscillatory forces, while the growth in the other 
two directions is slowed down. For surfactant concentra- 
tions p s < 0.2 anisotropic growth in x- and y direction is 
observed depicting the presence of tubular domains in the 
system. In future work it would be of interest to study 
the formation of tubes and their dependency on the shear 
rate, oscillation frequency and surfactant concentration 



in greater detail. Additionally, the study of asymmetric 
mixtures, where the concentrations of the two immiscible 
fluid species differs, remains unexplored. 
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